Detectable seismic consequences of the interaction of a primordial 

black hole with Earth 

Yang Luo^ , Shravan Hanasoge^'^ , Jeroen Tromp^'^* and Frans Pretorius^ 



(N 

o 

(N 



^ Department of Geosciences, Princeton University, Princeton, NJ 08544, USA 

^ Max-Planck-Institut fiir Sonnensystemforschung, Katlenburg-Lindau, Germany 

^ Program for Applied & Computational Mathematics, Princeton University, Princeton, NJ 08544, USA 

** Department of Physics, Princeton University, Princeton, NJ 08544, USA 

* e-mail: jtromp@princeton.edu 



\o 



O 
U 

6 

u 



> 

o 

00 

en 
rn 

o 



X 



Abstract 

Galaxies observed today are likely to have evolved from density perturbations in the early universe. Perturbations 
that exceeded some critical threshold are conjectured to have undergone gravitational collapse to form primordial black 
holes (PBHs) at a range of masses { jZel'Dovich & Novikov|[T967l|Haw"king||1971||Carr & Hawkingl[T?74| l. Such PBHs 
serve as candidates for cold dark matter jCarr et al. 2010) and their detection would shed light on conditions in the 
early universe. Here we propose a mechanism to search for transits of PBHs through/nearby Earth by studying the 
associated seismic waves. Using a spectral-element method ( Komatitsch & Tromp|p999:) , we simulate and visualize this 
seismic wave field in Earth's interior. We predict the emergence of two unique signatures, namely, a wave that would 
arrive almost simultaneously everywhere on Earth's free surface and the excitation of unusual spheroidal modes with a 
characteristic frequency-spacing in free oscillation spectra. These qualitative characteristics are unaffected by the speed 
or proximity of the PBH trajectory. The seismic energy deposited by a proximal M^^^ — 10^^ g PBH is comparable 
to a magnitude Afw ~ 4 earthquake. The non-seismic collateral damage due to the actual impact of such small PBHs 
with Earth would be negligible. Unfortunately, the expected collision rate is very low even if PBHs constituted all of 
dark matter, at ~ 10~^ yr^^, and since the rate scales as 1/A/^^^, fortunately encounters with larger, Earth-threatening 
PBHs are exceedingly unlikely. However, the rate at which non-colliding close encounters of PBHs could be detected by 
seismic activity alone is roughly two orders of magnitude larger — that is once every hundred thousand years — than the 
direct collision rate. 

1 Introduction 

The determination of the constituents of dark matter is an outstanding problem that continues to receive great attention 
2010[l. Constraints derived from the gravitational microlensing survey Experience de Recherche d'Objets Sombres 



(Feng 



C. Alcock et al. 1998 |P. Tisserand et al. 2007| l limit at most 8% of our galaxy's dark matter halo to be comprised of 
objects with mass M^^^ between 10'^'' g and 2.4 x 10"^^ g (for comparison. Earth's mass is approximately 6 x 10^^ g). 



< 10^^ g would emit 
2010|l. Based on these 



PBHs, which co uld constitute a significant fraction of these massive objects, are subject to Hawking radiation losses 
(Hawking 1974 1, a process by which they emit photons of a characteristic frequency ex 1/A/^^^. PBHs with masses 
]\/[fnti ^ lO^F^g would thus evaporate within the age of the universe, and those with M^^^ 
potentially detectable Gamma rays ( |Page & Hawking 1976; P. Sreekumar et al. 1998 Carr et al. 

criteria, PBHs in the mass range 10^^ g < M^^'^ < 10^^ g are unconstrained and the entirety of dark matter in our galaxy 
could in principle be constituted by them. One way to narrow these constraints is to search for special oscillations of the 
Sun (and potentially other stars) induced by the gravitational interaction with proximal and impacting PBHs. The impacts 
of PBHs of the size of asteroids (A/^^^ = 10^^ g) with the Sun would be detectable by existing solar observatori es, but 
the rate of occurrence of such events is estimated at a disappointingly low 10^® yr^^ ( |Kesden & Hanasoge 201 1 1. 

The prospect of black holes striking Earth excites considerable public interest. The PBHs we consider here are of 
masses small enough that they do not pose an existential threat to our planet, and because they are at the lower end of 
the allowable mass range, will cause negligible non-seismic effects. However, for masses the size of asteroids and larger, 
devastation similar in some manner to the impact of an asteroid may ensue. Consider, for example, the Tunguska impact 
and explosion in Siberia on June 30, 1908. No comet ash has ever been recovered around this region, leading to the 



speculation that it might have occurred due to an impact with a PBH (Jackson & Ryan 1973 1. Subsequently, this theory 



was dismissed because of (among other reasons) the lack of secondary signatures around the exit location of the supposed 



PBH ( Beasley & Tinsley 1974, Burns et al. 1976a I. Another example is the putative identification of supersonic waves in 



recorded seismograms, hypothesized to have been caused by the impact of a 'quark nugget' with Earth ( Anderson et al. 



[2003 ^ on Nov. 24, 1993. This claim was later falsified when it was shown that the signals were actually due to a hitherto 
unreported earthquake ( |Selby et al.| |2004) l. Quark nuggets (QNs), also termed strangelets, are a hypothetical phase of 
matter comprised of roughly equal parts of up, down and strange quarks that have lower binding energy than iron, and 
thus may be the "ultimate" stable state of matter ("Witten 19841. Because a QN is comprised of materials at nuclear 
density and can be of any mass up to a neutron star mass, it will have a very similar gravitational effect as a PBH when 
passing through Earth, and present a similar seismic signature. However, smaller QNs, of order M — \Q^^ g, could still 
produce observable explosions in the atmosphere, in contrast to PBHs of the same mass. This could be used to distinguish 
these two classes of compact object collisio ns with Earth. 

The speed of PBHs in unbound orbits (Lisanti & Spergel 201 l| l is on the order of 200-400 km s^^, indicating a 
maximum transit time through Earth of less than a minute. During its passage, a PBH would exert gravitational forces 
on the whole globe and generate seismic waves that propagate through Earth's interior. Due to the fast movement of 
the PBH compared to typical seismic wave speeds, such an event may be treated effectively as a supersonic source that 
generates Mach waves. The signatures of Mach waves are different from those of earthquake-induced seismic waves, 
thereby providing a means of detecting PBHs. Amplitudes of such waves have also been estimated ( [Khriplovich et al.] 
2008 1, albeit in a simplistic scenario where Earth is treated as a homogeneous fluid ball. 

Gravitational forces in the near-field of the PBH are very strong, leading to non-linear effects such as the cracking 
of rocks and heating. However, in this paper we focus on the far-field, treating the PBH as a moving Newtonian grav- 
itational source to the seismic wave equation ( [Kesden & Hanasoge[ |2011[ l, which we solve using the publicly available 
SPECFEM3D_GL0BE software package. An outline of the rest of the paper is as follows. In Section |2] we introduce 
the governing equations and describe the numerical solution method. In Sections [3] and [4] we discuss the results, in the 
latter section focusing on the unique spherical normal modes induced by a PBH. In Section l5] we provide some order- 
of-magnitude estimates of event rates and non-seismic effects. In Section l6] we suggest a strategy for detection, before 
concluding in SectionJT] We relegate to the appendix convergence tests and a few additional simulation results. 



2 Governing Equations and Numerical Implementation 
2.1 Governing Equations 

In seismology, the displacement field s of particles within Earth's volume Vl is governed by the seismic wave equation 

pa^s-v-T-f , (1) 

where p is the mass density, T the symmetric second-order stress tensor and f the source. Equation (fill is a linearized 
version of conservation of linear momentum, which proves to be sufficiently accurate for seismic wave propagation in 
Earth's interior 

In an elastic material, the stress tensor T is determined in terms of the strain via Hooke's law. 



T = c :£ 



(2) 



where c is the fourth-order elastic tensor with at most 21 independent components, describing the elastic properties of the 
material, and e the strain tensor 



e = 



Vs + (Vs) 



(3) 



In an isotropic elastic medium, the fourth-order tensor reduces to 2 independent parameters — the bulk modulus k and 
the shear modulus fi. Furthermore, the shear modulus /i vanishes in an acoustic medium (e.g., the outer core). Therefore, 
only the parameter k is relevant in fluids. It is important to note that because of this difference between solids (inner core, 
crust & mantle) and liquids (outer core), the governing equation takes sUghtly different forms, as will be further discussed 
later. Complications associated with attenuation are readily captured by an absorption band model ( jDahlen & Trompl 
1998) , and are accommodated in numerical simulations based on the introduction of memory variables (e.g., Komatitsch] 
& Trompl [T999l l. 

In global/regional seismology, the source term f represents earthquakes, whereas in exploration seismology, active 
sources, such as man-made explosions, are usually employed. To study the effects of PBHs, neither source type is 



relevant. Following Kesden & Hanasoge ( 2011)l, we consider a classical black hole, for which the source term f may be 



written as the gradient of the gravitational potential (p of the black hole: 

f= -pV(/. . (4) 

The potential (j) obeys Newton's law of gravitation 

_ GMPBh 

^^"^^ ^>~ ~ ||x-xPBH(i)|| ' (5) 

where G is the universal gravitational constant and A/^^^ the mass of the PBH. Note that the gravitational potential 
(hence the gravitational force) changes with time as the PBH is fast moving with a trajectory 

^PBH/.\ ^PBH ,,PBH. /^s, 

where Xq ^^ is its initial position (when we start the simulation) and v^^^ the velocity of the PBH. The energy loss of 
the PBH during its passage through Earth is so small that the velocity of the PBH remains constant. Consequently, the 
source term f takes the form 

f(x,^)=-GMPBH,(,) ^^^:5bh[^^)||3 ■ (7) 

Given its mass A/^^^ and trajectory (roll, the gravitational force that the PBH exerts on particles at any location within 
Earth may be calculated. 

On Earth's surface dVl, traction-free boundary conditions must be satisfied, 

n • T = (x e ar^, Vt) , (8) 

where n denotes the unit outward normal along the boundary, i.e., on the free surface dVt in this case. It is critical to 
note that Earth consists of several major shells: the crust and mantle that behave as solids on seismological time scales; 
the outer core that freely flows like a liquid; the inner core that is a dense solid comprised of mostly iron. Hence, the 
boundaries between these layers, namely, the Core-Mantle Boundary (CMB) and the Inner-Core Boundary (ICB), involve 
solid-fluid interactions. The boundary conditions on fluid-solid internal boundaries are continuity of normal displacement 
and traction, 

[n • s]+ = (x e S, Vt) , (9) 

[n • T]; =. (x e E, \/t) , (10) 

where [g]_ denotes the jump of any function g across an internal boundary E (either the CMB or ICB). Together with 
the initial conditions 

s(x,0) = O , ats(x,0)=O , (11) 

the system of governing equations presented in this section may now be solved numerically. 

2.2 Weak Form 

The wave equation (fill is in differential form, which is generally referred to as the strong formulation. Classical numerical 
methods, such as finite-difference and pseudo-spectral methods, usually involve a strong formulation. However, to account 
for solid-fluid boundary conditions properly, a weak formulation or integral form of the wave equation is preferable. The 
weak formulation is obtained by taking the dot product of equation ([TJ with an arbitrary test vector w, and integrating 
over Earth's volume Vt: 

(12) 



( (p 9^3 _ V . X) . w d^x = /" f • w d^x 
Jo, Jn 



2.2.1 Crust & Mantle 



The crust and mantle occupy the solid outer parts of Earth. Using the divergence theorem, equation ( [T2| in the crust and 
mantle may be simplified as 

/ pd'^s-vid^x+ Vw:Td3x-/ n • T • w d^x + / n • T • w d^x = / f • w d^x , (13) 

JQm -J^m JdQ. JcMB Jum 



where f^M denotes the volume of the crust and mantle, CMB the Core-Mantle Boundary and dfl Earth's free surface. 



The traction continuity condition on the CMB, i.e., equation ( 10 1, reduces to 



hT = -ph (x e CMB, Vi) , (14) 

where p is the pressure field in the outer core. 

Taking into consideration the free surface boundary condition (|8]l and traction continuity condition on the CMB ( [T4| , 
equation ( [T3] l may be rewritten as 



pdts-wd-^x+ Vw:Td-'x-/ p n • w d^x = / f • w d-'x . (15) 

It is the pressure field on the CMB that facilitates the exchange of information from the outer core to the mantle. 

2.2.2 Outer Core 

The outer core is a liquid in which the governing equation simplifies significantly. The di splacement field in the outer core 

^fluid ' ' 



may be expressed in terms of a scalar potential x as (Komatitsch & Tromp 



2002a I 



^fluid _ 1 



Vx . (16) 

P 



Substituting this definition into equation (fTJ, we may obtain a scalar governing equation for fluids 



K \p J K 

The weak form of this equation is 



^-92^-V-(-Vx) =--0 . (17) 



/ —dfx^^^+l -Vw-Vxd^x-/ — Vx-nd^x+/ — Vx ■ n d^x = - / — 0d^x 



(18) 



where JIq denotes the volume of the outer core and w is an arbitrary scalar test function. The continuity condition for the 
normal component of the displacement field on the CMB and ICB Q reduces to 

-Vx-n = n-s"°'''' (x e CMB or ICB, Vi) , (19) 



where s^°'"^ refers to the displacement on the solid side of the boundary. Now equation ( 18 i may be rewritten as 

/ — 9t^xd^x+/ -Vw-Vxd^x-/ wn-sd^x+/ w n • s d^x == - / — 0d^x . (20) 

Juo '^ "'Oo P JcMB JlCB Jno '^ 

It is the normal component of the displacement field at the CMB & ICB that passes information from the mantle & the 
inner core to the outer core, respectively. By definition, the pressure field in the outer core is 

p= -kV-s= -kV- Qvx) = -dh-cl> , (21) 



where in the last equality we have used equation (17 1 



2.2.3 Inner Core 

The weak form in the solid inner core is similar to the crust and mantle, namely 

p dtS • w d^x + / Vw : T d^x + / p n • w d^x = / f • w d^x , (22) 

Oi ^Oi JicB Jni 

where fti denotes the volume of the inner core. Equations ([TSj, pO] ) and (|22]) only hold under simple circumstances, be- 
cause they do not take into account rotation or self-gravitation. General but more complex formulations may be found in 
[Komatitsch & Tromp|(|2002a|b|l; such complications are accommodated by the spectral-element solver SPECFEM3D_GLOBE. 



2.3 Black Hole Singularity 

One of the main differences between earthquake and black hole simulations lies in the fact that we have to deal with 



different source terms in equations ( 15 1, (20i and (22 1, namely 



— w d X and 



f • w d X 



Om+^^I 



where and f are determined by dSl) and (JTl), respectively. Unlike earthquake simulations, where sources are associated 
with fault planes, here we encounter volumetric body forces. When a point within Earth is close to the location of the 
black hole, the gravitational force and potential become singular To avoid such unphysical singularities, we implement 
the source term as 



,PBH 



f(x,t) 



GA'r''''p{x 



G M'-'''' p{i^) 



it) 



^PBH 



(tW 



,PBH 



it) 



-' '■crit 



cP^"(t)||>i?crit , 



(23) 



GM 



PBH 



(x,i) 



,PBH 



it)\\ 



3GMPBH gmPbh 



PBH 



2R. 



'crit 



2 /?3 

^ ^crit 



it)\\' 



|x-xPBH(i)||>^^^.^ 
|x-xPBH(i)||<^, 



'crit 



(24) 



where J?crit is a chosen critical length scale within which the amplitude of the gravitational force decreases linearly. It 
may be shown mathematically that and f capture all features of the original source terms and f , as long as the critical 
length is less than the size of the numerical grid. In our simulations, the size of an element is ^ 30 km, while i?ciit is 
around ^ 6 m, i.e., five thousand times smaller than the element size. Convergence tests using different i?ciit illustrate 
that further decreasing i?ciit does not change the simulation results (Appendix [A]). 

The last, but equally important issue is the selection of the time step for the numerical simulations. As mentioned 
previously, the speed of the black hole is typically a hundred times larger than the seismic wave speeds, and it takes only 
a few minutes for the PBH to traverse Earth. Therefore, we need a much smaller time step than usual to resolve the 
trajectory of the black hole. Fig.|9]in Appendix |A] explains why this is important. In practice, we choose a time step that 
is one thousandth of the original time step when the PBH is close to Earth. 



3 Results of Spectral-Element Simulations 

The set-up of our simulations is illustrated in Fig.[T] Different trajectories as well as two distinct black hole s peeds, namely 



Komatitsch & 



200 km s~^ and 400 km s~^, are considered. The simulations are carried out on a cubed-sphere mesh ( 
Tromp 2002a|b I comprising 8.2 million elements. Given the wave-speed structure of Earth, the shortest period the mesh 
can resolve is around 15 s. We set the black hole mass Af ^^^ — 10^^ g, noting that amplitudes of the resultant seismic 
waves scale linearly with A/^^^. 

An instructive case is that of a PBH whose trajectory coincides with the center of Earth, i.e., li = in Fig.[T] We use 
the 1-D Preliminary Reference Earth Model ( [Dziewonski & Anderson[[1981[ ), aligning the PBH trajectory with Earth's 
rotation axis. Four snapshots from this simulation are shown in Fig. |2] illustrating generation, propagation, transmission 
and reflection of PBH-induced seismic waves. Seismograms from this simulation are displayed in Fig. l3] Due to sym- 
metry, East-West ground motions are very weak compared to the other two components, and are thus not shown. The 
North-South component shows not only distinctive arrivals due to cylindrical waves from the supersonic motion of the 
PBH, but also classical Rayleigh surface waves and reflected/transmitted core phases. The vertical component, besides 
containing signatures similar to the North-South component, registers unique waves that arrive at almost the same time 
everywhere on Earth's surface, distinct from those arising due to earthquakes. These waves are generated at the Core- 
Mantle Boundary (CMB) and Inner-Core Boundary (ICB) due to free-slip interactions between the solid mantle & inner 
core and the liquid outer core. 

We study four additional PBH trajectories, two of which pass through Earth's interior (d — 0.4 i?® and d — 0.8 i?®) 
while the other two pass nearby (d = 1.2 R^ and d = 1.6 i?©). A faster PBH traveling at a speed of 400 km s~^ is 
also considered. To avoid redundancy, only d = 0.8 i?© for the faster PBH is shown in Fig. HI (more images may be 
found in AppendixlBll. On the North-South component, clearly visible are direct arrivals due to the supersonic motion of 




Figure 1: Set-up for the numerical simulations. The thick black line denotes the PBH trajectory, where the distance 
between the trajectory and the parallel diameter of the Earth is denoted by d. For example, when the PBH passes through 
the center of Earth, d = 0; when the PBH misses the Earth, d > i?®, where i?® — 6371 km is the radius of Earth. 
Hypothetical seismic stations denoted by triangles are deployed around the globe, but only the most significant ones 
(yellow triangles) will be shown later. 
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Figure 2: Snapshots at four different instantaneous times. Shown are particle motions along the a:-direction in the x-z 
plane through Earth's center Red represents motion to the left, i.e., the negative x direction, whereas blue indicates that 
particles are moving to the right, a Generation of cylindrical and spherical waves. The PBH has entered Earth from 
the top and is about to exit from the opposite side. Some wave groups are labeled using commonly used seismological 
nomenclature, and may be viewed together with those shown in Fig. 3. b Propagation of cylindrical and spherical waves. 
Spherical waves from the CMB reach Earth's surface everywhere at the same time, leaving the linear features we see 
in Fig. 3. Because of the non-homogeneous seismic wave speeds, cylindrical waves are bent and gradually lose their 
cylindrical shape, c Spherical waves bounce back from the surface, whereas cylindrical waves continue to arrive at 
stations from high latitudes to the equator, d Spherical waves from the ICB reach the surface and cylindrical waves 
generated in the outer core & inner core (KP/IKP) are about to hit the equator. It is important to understand that these 
snapshots are for illustrative purposes only. A large critical length scale is used to eliminate high-frequency numerical 
"noise" and hence the amplitudes of cylindrical waves are greatly underestimated compared to spherical waves. 
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Figure 3: Simulation results for v = 200 km s^^ and d = in Fig. 1. East- West components are omitted since no 
signals can be observed due to symmetry. Snapshots in Fig. 2 may help to better understand phase labels, a North-South 
components of ground motion. The earliest arrival (solid blue line) is often called the Mach wave, similar to an earthquake 
P phase. The solid red line represents seismic waves reflected at the CMB, corresponding to an earthquake PcP phase. 
The solid green line highlights the PKP phase traveling through the outer core. Phases KP and IKP are represented by 
the dashed blue line, generated either in the outer or inner core, respectively. The KKP phase, denoted by the dashed red 
line, is reflected once at the CMB before entering the mantle. The phase denoted by the dashed green line has a wave 
speed of ^ 3.8 km s^^, i.e., the Rayleigh surface wave, b Vertical components of ground motion. Besides imprints of 
the same phases seen on the North-South components, phases emerge that arrive at all stations simultaneously. These 
involve spherical boundary waves generated at the CMB & ICB, labeled by dashed blue and red lines, respectively. Time 
intervals between these phases correspond to reverberations. For example, the time interval between the first and second 
blue dashed lines corresponds to waves traveling across the outer and inner cores through the Earth's center, while the 
time interval between the first and third blue dashed lines corresponds to waves traveling across the entire planet. 
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Figure 4: Simulation results for d = 0.8 i?® and v = 400 km s~ . a North-South component ground motion, b Vertical- 
component ground motion. Both cylindrical and spherical waves can be seen clearly. Focusing points of the Rayleigh 
surface waves move to latitudes of ±40°, reflecting changes in the entrance and exit locations of the PBH. 



the source as wefl as Rayleigh surface waves, which behave in the same manner as those shown in Fig. [3] Distinctions 
from the simple example also emerge. For example, the entrance and exit points of the PBH change from the poles to 
latitudes around ±40°, as indicated by the focus of the surface waves. Waves arriving after the surface waves are coherent 
but more difficult to interpret. The vertical component is largely unchanged from the simple example. Besides other 
signatures found on the North-South component, linear features due to the spherical waves exist as well. For cases where 
the PBH does not hit our planet, we see no coherent signals on the horizontal components. But spherical waves are still 
visible on the vertical component, since they are gravitationally induced at the CMB and ICB. Simulations in 3-D Earth 
model S40RTS ( [Ritsema et al.|[20TT| l preserve these key features, as illustrated in Appendix [B] 

For all cases, ground motions are on the order of 10^^ m. For comparison, ground motions due to the Bolivia 
il/w = 8.1 earthquake on June 9, 1994 are on the order of lO^'^ m. From an energy perspective, seismic energy injected 
by the 10^^ g PBH is six orders in magnitude lower than that of the Bolivia earthquake, corresponding to an earthquake 
of magnitude M^, = 4. 
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Figure 5: Stacked spectra of vertical-component seismograms at all stations. Long time series of 20 hours are used to 
calculate the spectra, instead of the 50-minute seismograms shown in Fig. 3. Labeled normal modes are from a mode 
calculation in 1-D model PREM, where red lines denote modes that are not commonly observed, such as 25*2 ■ Differences 
between the peaks and the labeled modes are ^ 0.02 mHz, around the sampling frequency resolved by a 20-hour long 
record. Besides these unusual modes, the evenly-spaced peaks define another distinct PBH feature. The behavior of these 
modes is similar to radial modes, „S'o, and their frequency spacing corresponds to the compressional wave travel time at 
an epicentral distance of 180° (PKIKP). These modes correspond to the spherical boundary waves we see in Fig. 3. 



4 Normal-Mode Excitation 

Another distinctive feature of PBH-induced seismic waves is the appearance of unusual spheroidal normal modes (which 
are not excited by earthquakes) in spectra of the time series, as shown in Fig. l5] 



The significance of the power spectra arises from normal-mode analysis (Dahlen & Tromp 19981. Any motions on 



Earth can be decomposed into a set of periodic oscillations — Earth's free oscillations or normal modes: 

s(x,i) = Rc^ (iwfe)-isfe(x) /" Ck{t')eyi^[iuJk{t-t')]At' , 
k "'-°° 

where s^ denotes the fcth mode, ujk its eigenfrequency and Ck its excitation. 

4.1 Earthquake Excitation 

For an earthquake, the excitation term Ck is controlled by the location of the earthquake and its mechanism: 



(25) 



Cfe(i)= / m(x,t):£t.(x)d2x 



(26) 



where E* represents the fault plane at time t, m the moment density tensor and Sk the strain tensor associated with the 
fcth mode 



1 r. 



£k 



(27) 



When the point source approximation applies, as is often the case in global seismology, the excitation term reduces to 

Cfc(t)«M(t):4(x,) , (28) 
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where the strain tensor Sk at the location of the point source x^ interacts with the moment tensor 

M(t) = / m(x,t)d2x . 



(29) 



Equation ( 28 1 determines not only which modes can be excited by an earthquake, but also how strongly excited these 



modes are. It is critical to note that earthquakes only occur in the crust & upper mantle. Therefore, modes that are not 
sensitive to the shallow parts of Earth cannot be excited by natural earthquakes; examples include the Slichter mode (i5i) 
and Stoneley modes trapped at the CMB & ICB. 

4.2 PBH Excitation 

A PBH triggers normal modes in a different way than earthquakes. The excitation follows the form of a distributed force 



where f satisfies 



Ck{t)= [ f(x,i)-s^(x)d3x , 
Jn 

f(x,i)^GMPB"p(x) V^ 



,PBH 



Wll 



(30) 



(31) 



Substituting equation ( 31 1 into ( 30 1, we obtain 



Cfe(t) = / GMPB«p(x)s:(x).V 
Jn 



-. oo 
r' 2-^ 



An 



r — V 2/ + 1 
(=0 ^ 



r 



I 



-) Y. yun{e,4>)YL{e',c^') 



G iVr^" p(x) s^(x) ■[vdr + r-'Vi) 



rn^ — l 

An 



1=0 ^ 



d'y 



^ OO 



;) E ^""( 



y, 



Im.V' :' 



rn^ — l 



d3> 



(32) 



fiml. 



where r = 1 1 x 1 1 , r' = 1 1 x^^^ 1 1 , V i = 9 dg + cf) (sin 9)^^d^ denotes the surface gradient on the unit sphere, and Yj, 
are spherical harmonics (Edmonds 1960 1. Note that this expansion is valid only for r' > r. When the PBH is within 
Earth, r may be larger than r' for some points x. In that case, we need to interchange r and r', leading to a more complex 
result than the one considered here. 

Mode Sf; may be a spheroidal mode '^Sj. or a toroidal mode ^Sj,: 



^Sfe(x) = 

^s,(x) = 



«c//(r)r,™(0»f + 

-1 



1 



^/W+T) 



yi(ITT) 
nWi{r) (f X Vi)y^„ 



nViir) ViFi, 



(33) 
(34) 



where I, m and n are the degree, order and overtone number of the mode, respectively; „C//, „V; and „Wi are the 
corresponding radial eigenfunctions. Due to the relations 



f-Vi = 

f-(fxVi) = 

Vi-(fxVi) = 



(35) 
(36) 
(37) 



it is obvious that no toroidal modes can be generated by a PBH, because the excitation term vanishes. For spheroidal 
modes, equation ([32| reduces to 



Cfc(i) 

nAlit) 



An a A/P^H Z"^® r , 

(2,.i,(.mr i ^^') '■'^'[' Mr), .MTY) Miry 



dr 



(38) 
(39) 



where we have used the orthogonality of spherical harmonics. Note that the PBH location x^^^ varies with time, and 
therefore r', 9' and 0' are also functions of time. Equations ( 38 1 and ( 39 1 differ from equation ( 28 i in that the eigenfunc- 



tions everywhere along Earth's radius are involved. The implication is that even modes whose sensitivities are restricted 
to the deeper parts of Earth, i.e., modes that cannot be triggered by earthquakes, may be excited by a PBH. 
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Substituting equations ([38]) and ([39| into equation ([25]l, we obtain 



s(x, t) — ReTJ 



iWfc) 



-1 s 



Sfe(x) 



.^i(i') >^;:n ( e'{t'),<j>'{t') ) exp[zc^fc(t - t')] dt' 



(40) 



Note again that the summation over k actually denotes triple summations over I, m and n. In order to evaluate equation 
(pOb efficiently, the summation over m is carried out analytically using the spherical-harmonic addition theorem 



^ Y;^i0',cj,')Yiraie,c 



-I 



21 + 1 

47r 



PiicosO) , 



(41) 



where cos0 = cos cos 6*' + sin6'sin6''cos((/) — (f)')- Finally, equation (40 1 reduces to 

-IRoV (i^LoA-^ ^UAr\r f 
in 



s(x,i) 



l,n 



Ai{t')Pi{Q{t')) cxpH „wi t') dt' 



1 



nVi{r) / nMt') ViPi (e(i')) exp(-i ^uji f) At' 



exp(i n^i t) 



(42) 



.JW+l) 
Equation ([42]i is difficult to implement for a moving PBH, however, it is straightforward to consider a stationary and 



transient PBH effect. In this case, instead of equation ( 3 1 1, the force term is 

PBH ■ - — 1 



f(x,i) = C?Ar""p(x) V 



^PBHl 



5'{t) , 



(43) 



where 5'{t) — rather than 5{t) — is used to avoid inducing a permanent displacement. In this simpler case, where r', 9' 
and (j)' are time-independent, equation ( 42 1 takes the following form 



«(x,i)=(-^)Re^ ^Ai 



^Uiir) Pi (cos Q) r 



1 



^/WTT) 



,Vi{r)'^iPi{cosQ) 



exp(i riUJi t) . (44) 



We shall use this configuration solely for benchmarking purposes, testing the generation of spherical waves by comparing 
the results of a normal mode calculation with those obtained based on a SEM calculation. For earthquake sources and 
impacts, the SEM has already been benchmarked extensively against normal modes (e.g., |Komatitsch & Tromp| |2002at 
Meschede et al. 2011 1. Figure [6] shows a comparison between a normal-mode calculation and a spectral-element simu- 
lation for two stationary and transient PBHs, which are symmetrically placed with respect to Earth's center. Using two 
PBHs instead of one ensures that Earth's center remains stationary. 



5 Estimated Event Rates and Non-seismic Effects 

In this section we provide some order-of-magnitude estimates of event rates and non-seismic effects of a black hole (BH) 
or quark nugget (QN) moving through the earth. The radius of a BH is 



TBH ~ 1-5 X 10^'^mM 



(45) 



where M = M/10^^ g is its mass in units of 10^^ g. Assuming a QN is a uniform sphere with density pqn 
lO^^g/m'^po, with po of order unity (Witten 1984 1, it has a radius of 



TQN « 6 X 10 m 



M 



1/3 



(46) 



Restricting the discussion to masses less than ^ 10~^M(^, one can see that though the QN is significantly larger than the 
BH, as gravitational sources both can still be considered point-particles to good approximation. Assuming that a fraction 
X of the estimated dark matter density in the vicinity of the solar system — see for example Lisanti & Spergel ( 2011[ l — 
is in the form of compact objects of mass M, the expected event rate is 



iV « 4 X IQ-Vr-^ ^^ 
M 



(47) 
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2000 
Time (s) 



Figure 6: Benchmark between a normal-mode calculation and a spectral-element simulation for two stationary, transient 
PBHs placed symmetrically about Earth's center. This benchmark confirms the accurate implementation of the volumetric 
PBH source in the spectral-element code. 



where v = t;/200 km/s is the velocity in units of 200 km/s, and d = bainx/R® is the maximum distance, in units of i?^, 
out to which the passage of a compact object with impact parameter femax will create the kind of impulsive gravitational 
response that gives rise to the unique seismic signatures discussed here. A rough estimate suggests d w 10 — 20. These 
events are thus quite rare for masses that could produce observable seismic signatures. However, as we will show below, 
on the lower mass end of BHs that are seismically observable, the corresponding non-seismic signatures will be negligible, 
and seismology may offer a unique way to observe them. This is not the case for QNs, which, because of their larger radii, 
could produce strong atmospheric explosions. 

The Bondi radius Rb for accretion of material onto a compact object is 



Rb ~3x lO^Vi 



M 



(48) 



For the masses considered here, this is negligibly small for BHs, and is inside the radius of a QN, thus for QNs the relevant 



radius is simply tqn. Following Bums et al. U976b[ ), the radius Rm out to which the gravitational impulse of the passing 
compact object could melt surface rock is 

R,n~5x 10" V — , (49) 



and following |Labun et al.| ( |2011[ ), the transverse distance L out to which rock will be shattered due to the gravitational 
tidal field is 



10"^m M 



(50) 



Again, both Rm and L are within the radius of the QN, and neghgible for smaller mass BHs. Applying the results 



of Ruderman & Spiegel ( 1971 1, we may estimate the energy that will be deposited into the atmosphere by the passage of 
the compact object as it passes through. For BHs the main effect is due to the impulsive gravitational acceleration of air 
molecules, giving 

AE^%^^0.3 3 ( ^j . (51) 

To create an explosion such as the Tunguska event of ^ 10^^ — 10^^ J, for example, would require a BH of mass 
10^ — 10^ M. QNs cause a similar gravitational impulse, however their larger surfaces experience a drag that can deposit 
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significantly larger amounts of energy compared to BHs for the smaller mass objects: 



AE^^' « 4 X 10^° J vM — J . (52) 

Again, to create Tunguska-like explosions requires massive QNs, however, in contrast to BHs, smaller QNs can still 
produce sizable explosions, e.g., a M — 1 QN imparts the energy equivalent of ^ 10 tons of TNT. 



The event rates above are exceedingly low for massive (Af ^ 1) compact objects, though as discussed in Labun 
|et al.|pOTT) , direct impacts may leave features distinguishable from asteroid impacts, and hence the effective baseline 
for observation can be extended by looking at geological features on Earth, Moon, Mercury and Mars. Also, as BH and 
QN detectors the gas giants have larger collecting areas, and may be useful in constraining compact object dark matter 
candidates in regimes not readily accessible by other observations. 

6 Recipe for Detection 

Based on numerical simulations, we find that PBH transits through Earth generate unique and significant seismic signals. 
Our theoretical results may be used to develop a strategy for detecting PBHs. Unlike earlier suggestions (Anderson et alj 
2003|l, we propose to use spherical rather than cylindrical waves as the key detection mechanism. There are several 



advantages to this approach. First, detection is more reliable because the simultaneous global arrivals induced by slip at 
the CMB and ICB are unique and cannot be produced by earthquakes. In contrast, cylindrical-wave arrivals may resemble 
certain earthquakes, especially when the number of available stations is small. Such a misinterpretation, involving an 
unreported earthquake, occurred previously (Selby et al. 2004)). Second, in order to detect cylindrical waves, travel time 



measurements are involved, which presents a by-no-means trivial task. When the trajectory is unknown, it is difficult to 
align seismic stations, which further increases the challenge of picking the relevant arrivals. In contrast, the linear travel- 
time feature associated with simultaneous global arrivals can be simply revealed by comparing seismograms side by side 
in a so-called record section. Because the signals have the same arrival time, how we align the related seismograms is 
less critical. Third, due to reverberations of the CMB- and ICB-induced spherical waves, unique free-oscillation spectra 
are expected. These spectra involve a frequency-spacing between modes excited by a PBH governed by the pole-to- 
pole compressional-wave travel time, and the excitation of modes not affected by earthquakes, which are restricted to 
Earth's crust and upper mantle. Stacking of global seismic spectra may be required to bring out these free-oscillation 
characteristics. Finally, the detector we propose applies equally to scenarios where the PBH does not collide with Earth, 
but passes sufficiently close to generate detectable spherical waves. Of course cylindrical waves are still important, 
because fitting the arrival times of Mach waves provides constraints on the PBH trajectory which cannot be obtained by 
considering spherical waves only. 

Seismograms recorded by the global seismographic network are routinely monitored to detect landslides and glacial 
earthquakes ( jEkstrom et al.|[2003] l, and it may be feasible to include detection criteria for PBHs in this monitoring process. 

7 Conclusions 

We study the effects of PBHs passing through or nearby Earth using numerical simulations based on a spectral-element 
method. The fast motion of the PBH causes rapidly varying classical gravitational forces that are carefully numerically 
accommodated. We also satisfactorily resolve numerical issues due to the singularity in the gravitational near-field of 
the PBH. we predict three characteristic features. The first feature involves cylindrical waves due to a supersonic source, 
generating Mach waves which are readily understood and have been studied previously ([Anderson et al. 2003 1. The 



second feature involves spherical waves generated by free-slip at solid-fluid boundaries, namely the core-mantle and 
inner-core boundaries. These spherical waves present a unique signature at global seismographic stations that cannot 
be produced by earthquakes. Third, free-oscillation spectra are predicted to exhibit evenly-spaced peaks, with a spacing 
determined by the pole-to-pole compressional-wave travel time. Modes not excited by earthquakes, such as 2>5'2, should 
be detectable in spectra induced by PBHs. Probably the most likely scenario is the passage of a PBH nearby Earth. In 
such a near-miss scenario, the cylindrical-wave features vanish, but the spherical waves remain. Therefore, we may still 
use global seismograms and their spectra to detect a PBH, although it would be difficulty to determine its trajectory. We 
also note that spherical waves could in principle be excited by asteroids along trajectories proximal to Earth, potentially 
containing new seismic information about the interior 

In this article we only consider gravitational interactions between a PBH and Earth. Potential PBH impact phenomena, 
which are still debated, are not considered. Related seismic effects may be readily taken into consideration, as exemplified 



by the meteorite impact study of Meschede et al. ( 201Ij ). Such effects are of course irrelevant for near misses. 
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A Convergence Tests 



As discussed in Section 2.3 we introduce a critical length scale i?ciit in order to eliminate unphysical singularities. This 
length scale determines the maximum force we allow in the numerical simulations, as well as the behavior of the source 
time function. On one hand, amplitudes of simulated seismograms will be too small if i?crit is too large — the forces 
the PBH exerts on Earth would be severely underestimated. On the other hand, the smaller i?crit, the more the effective 
source time function looks like a delta function, which contains high-frequency information that cannot be resolved by our 
finite-size elements. The second issue is less problematic, since high-frequency noise can be safely removed by applying 
low-pass or band-pass filters in post processing. For these reasons, i?crit should be chosen as small as possible, as long as 
the maximum force does not exceed the limits of computer representation. However, we will show that when the critical 
length scale is far less than the size of the elements, differences in the weak form implementation vanish. 
Fig. |7] illustrates an element which contains the PBH at an arbitrary moment. 



Q. 




Figure 7: Illustration of how to avoid unphysical, singular gravita- 
tional forces from the PBH. The big square fte denotes the volume 
of an element and the circle represents a tiny sphere ilpBH around 
the PBH defined by the radius i?crit- Note that the cartoon is not 
drawn to scale — the size of the sphere is exaggerated for clear 
viewing. In fact, i?crit should be much smaller than the size of the 
element Ah. The net gravitational force within ilpBH vanishes, 
as long as i?ciit ^ ^h. Therefore, we can avoid singularities by 



Ah 



using the alternative force (23 i. The potential (24i is then derived 
from the corresponding force. 



The source term involves a volume integral of the gravitational force over element fie, that is 



f • w d'^x 



f • w d''x + / ( f - f ) • w d^x 



f • w d'^x 



G Af ^«" p(x 



x-x^i^"(t) x-x^«W(t) 



||x-xPBH(i)||3 ^3^.^^ 



f-wd^x+Z GAfPBHp(^PBH(^)) 

f • w d^x + G A/P«« p (xPBH(^)) ^ (xPBH(^)) 



wfx) d'^x 



x-x^i^"(t) x-xi^«W(t) 



(xPBH(^)) ^3^ 



X-xPBH(i)||3 ^3^.^ 

X-xPBH(^) x-xPBH(^) 

||x-xPBH(i)||3+^:^^ 



d^x 



where we have used the fact that when i?crit ^ Ah, flpBH can be treated as a single point compared to element fi^, i.e., 
p(x) w p {x^^^{t)) and w(x) « w (x^bh^^)) foj- ^ny x e flpBH at time t. 



Define r = x — x (f), and hence r 



X — X (i) 1 1 . Then the volume integral reduces to 



f • w d X - 



GAfPBH^^^PBH(^)) w(xPBH(;)) 






d^X 



(53) 



because the integral over fipBH is non-singular and simply vanishes. This can be easily understood from a physical 
perspective: within a small spherical volume around the PBH, the net gravitational force from the PBH is zero because 
gravity is isotropic in the tiny sphere. Equation p3| ) is important, because it provides a way to calculate the source term 
involving a singular force f using a non-singular force f , which is easily implemented in the spectral-element method. 
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We design an experiment to test this result. Three complementary simulations are carried out, each with a different 
i?crit: the first one with i?crit = lO^^i?©, the second one with Rait — lO^^i?© and the last one with i?crit = lO^^i?©, 
where i?® = 6371 km is the radius of Earth. One seismogram from each of the three simulations is shown in Fig. l8] 
The three simulations give the same results, thereby demonstrating that the critical length scale i?crit — lO^^i?© used in 
our simulations is sufficiently small for our theoretical derivation to be valid. Using f instead of f does not only removes 
singularities, but also gives reliable seismograms. 



xlO 



Convergence Tests 
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Figure 8: Convergence tests using three different i?crit- The perfect match shows that the critical length scale we used 
for our simulations is sufficiently small. Further decreasing i?ciit does not change the simulation results, as expected from 
equation ([53]l. 



Another issue involves capturing the widely differing time scales associated with the motion of the PBH and the 
periods of the induced seismic waves. Fig.l9]explains why this is important. 
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PBH 








Figure 9: Illustration of the significance of the choice of time step At. Two situations are shown, where the only 
difference is that At for the top row is twice that of the bottom row. The straight lines represent the trajectories of the 
PBH, whereas the crosses denote the PBH location at incremental time steps. The blue circles highlight the effective 
region of the PBH, i.e., the gravitational forces from the PBH outside these regions may be neglected because of the l/r^ 
fall-off of the gravitational field. Note that the radius of the effective region is not identical to the critical length scale 
^ciit- When At is chosen too large (top), there are points along the trajectory that are not sampled by the time steps. 
Therefore, the effect of the PBH is greatly underestimated. When At decreases (bottom), we cover the entire trajectory, 
which fully accounts for the gravitational effects of the PBH. 
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B More Examples 
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Figure 10: Seismograms from a simulation in which d = 0.8 i?© and v — 200 km s^^. Direct arrivals and Rayleigh 
surfaces waves are observed on both horizontal and vertical components, whereas phases that share the same travel time 
at all stations (straight vertical features in the figure on the right) only emerge on the vertical component. The focusing 
points of the surface waves are determined by the entrance and exit points of the PBH, while the different curvatures of 
first arrivals at high and low latitudes are an expected feature of Mach waves. 
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Figure 11: Seismograms from a simulation in which d — 1.6 i?® and v — 200 km s~^. The key observation for this 
scenario is that no coherent signals are observed on the North-South component, because the PBH does not transit Earth. 
When there is no direct contact between the PBH and Earth, the dominant mechanism changes from a supersonic source 
to differential gravitational forces at different locations in Earth's interior, a result which is more difficult to imagine and 
interpret. Nevertheless, not only does this result validate our explanations for the cylindrical- and spherical-wave features, 
it also provides a new way to detect a PBH, even when it passes nearby Earth. 
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Figure 12: Simulation result for 3-D Earth model S40RTS with d = 1.6 i?® and v — 400 km s ^. The characteristic 
signature from the spherical waves remains the same. 
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